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Plasma  Modeling" 

PI:  John  Loverich,  Tech-X  Corporation 
University  Collaborator:  Uri  Shumlak,  University  of  Washington 

Overview 

During  the  Phase  I  project  we  were  tasked  with  investigating  schemes  for  multi¬ 
fluid  plasmas  including  high  order  methods  on  unstructured  grids,  electric  field 
divergence  cleaning,  development  of  multi-fluid  benchmarks  and  the  development 
of  an  Oracle  to  tell  the  user  the  regimes  of  validity  of  various  plasma  models  in  a 
particular  simulation.  The  phase  I  project  was  quite  successful  and  details  of  our 
progress  are  described  below. 

Objectives: 

The  overall  goal  of  this  proposed  project  (Phase  I  and  II]  is  to  develop  a  high-fidelity,  efficient, 
accurate  and  easy-to-use  software  package  for  simulation  of  a  broad  variety  of  plasma  physics 
problems  relevant  to  Air-Force  applications.  Phase  I  was  meant  to  make  progress  in  this  direction. 

Status  of  Effort: 

Phase  I  is  finished 

Accomplishments/New  Findings: 

1)  Implemented  unstructured  grid  algorithms  (DG  and  finite  volume) 

2)  Added  improved  implicit  schemes  for  multi-fluid  modeling 

3)  Tested  new  divergence  preservation  approach 


Publications: 

As  part  of  a  general  description  of  USim  and  AIAA  paper  was  presented 
"Nautilus:  A  Tool  For  Modeling  Fluid  Plasmas",  lohn  Loverich.  Sean  C.  D  Zhou.  Kris 
Beckwith.  Madhusudhan  Kundrapu.  Mike  Loh.  Sudhakar  Mahalingam.  Peter 
Stoltz.  Ammar  Hakim.  51st  AIAA  Aerospace  Sciences  Meeting  including  the  New 
Horizons  Forum  and  Aerospace  Exposition.  2013, 10.2514/6.2013-1185 

This  paper  gave  a  general  overview  of  USim  [formerly  called  Nautilus)  in  its  current 
state  including  two-fluid  and  DG  algorithms  developed  during  this  Phase  I  project. 


Interactions /T  ransitions: 


Many  of  the  two-fluid  capabilities  developed  during  this  Phase  I  project  will  make  it 
to  the  first  release  of  USim  and  thus  be  available  to  the  public. 

New  Discoveries: 

No  new  discoveries 

Honors/A  wards: 

No  awards 

Task  1:  Create  Oracle  to  Determine  validity  of  models  for 
particular  applications 


Figure  1  shows  the  GEM  challenge  magnetic  reconnection  problem  with  the  Oracle 
applied.  Colors  indicate  the  regime  of  validity  of  each  of  the  plasma  models,  MHD, 
Hall  MHD,  two-fluid  and  then  kinetic  models.  It’s  interesting  to  note  how  a  problem 
that  can  essentially  be  described  by  two-fluid  or  Hall  MHD  initially  evolves  to 
something  where  only  kinetic  model  is  strictly  valid.  The  Oracle  has  been 
implemented  in  WarpX  and  will  be  implemented  in  USim  during  the  phase  II  project. 


Figure  1:  Oracle  guides  the  selection  of  appropriate  physics  model  which  are  nested 
sets,  i.e.  Kinetic  Contains  Two-Fluid  Contains  Hall-MHD  Contains  MHD.  Results  applied 
to  the  GEM  challenge  magnetic  reconnection  problem  with  nominal  values  from 
Earth’s  magnetosphere  identify  the  simplest,  applicable  model  for  the  initial  and  final 
simulation  conditions. 


Task  2:  Create  common  geometry  and  datastructure  library  for 
use  in  both  Nautilus  and  WARPX 


Tech-X  and  UW  are  currently  taking  different  approaches  on  the  unstructured  Grid. 
Both  organizations  have  finished  initial  unstructured  grid  implementations  and 


algorithms.  University  of  Washington  is  implementing  their  grid  for  GPU  and  MIC 
architectures  in  mind. 

This  unstructured  grid  architecture  in  USim  and  WarpX  allow  users  to  generate 
custom  meshes  with  increased  resolution  in  areas  where  small  scale  effects  such  as 
turbulence  is  expected,  or  meshes  that  need  to  conform  to  arbitrary  geometries.  The 
infrastructure  behind  mesh  generation,  conversion,  data  manipulation  and  domain 
decomposition  has  been  largely  completed  within  both  codes  for  ID,  2D  and  3D 
meshes.  Unstructured  meshes  are  currently  generated  from  the  mesh  generation 
suite  CUBIT  in  WarpX/M  and  using  Gmsh  in  USim  [addition  of  CUBIT  meshes  is 
occurring  as  part  of  a  separate  project  at  Tech-X).  In  WarpX/M  Python  scripts  are 
used  to  convert  the  mesh  files  and  generate  initial  conditions.  The  unstructured 
architecture  then  imports  the  mesh  and  distributes  the  workload  across  a  compute 
cluster.  Within  each  cluster  node  the  mesh  is  further  subdivided  into  patches  spread 
out  among  compute  devices,  such  as  CPUs  and  GPU's,  for  processing.  The  current 
task  is  to  utilize  this  infrastructure  in  the  implementation  onto  the  finite  volume  and 
finite  element  methods.  Figure  2  shows  a  flow  diagram  for  the  unstructured  mesh 
infdrasturcture  in  WarpX/M. 


Figure  2:  Process  diagram  of  the  unstructured  mesh  infrastructure  surrounding 
WARPX/M.  The  process  starts  by  designing  a  mesh  within  a  mesh  generation  suite  (Cubit). 
The  mesh  is  then  converted  and  read  into  WARPX/M  in  which  a  solution  is  generated. 
Solution  files  are  then  linked  to  the  mesh  file  through  an  xdmf  interface  file  and  analyzed  in 
an  unstructured  data  visualizer  such  Paraview  or  Visit. 

Task  3:  Evaluate  schemes  for  efficient  robust  solution  of  multi¬ 
fluid  equations 

There  are  3  major  problems  with  existing  Two-Fluid  algorithms  used  by  Tech-X  in 
Nautilus  and  the  UW  in  Warp-X.  (1)  a  good  divergence  preserving  approach,  (2) 
computational  restrictions  based  on  plasma  frequency  and  (3)  computational  restrictions 
based  on  the  light  speed. 

During  this  project  we  have  implemented  a  robust  solution  have  resolved  problem  (2)  by 
using  a  semi-implicit  two-fluid  (and  10-moment  2  fluid)  algorithm  designed  by  Harish 
Kumar  “Entropy  Stable  Numerical  Schemes  for  Two-Fluid  Plasma  Equations”  2012. 
Harish’s  work  is  similar  to  previous  work  performed  at  the  University  of  Washington 
however  it  uses  a  semi-implicit  Runge-Kutta  approach  that  preserves  high  order  time 
accuracy  while  being  an  unsplit  scheme.  This  allows  us  to  resolve  issues  that  we’ve 


observed  using  operator  splitting  and  lower  order  semi-implicit  integration.  Furthermore, 
the  approach  can  be  used  with  Discontinuous  Galerkin  since  it  does  not  require  operator 
splitting.  The  approach  works  for  both  the  two-fluid  and  ten-moment  plasma  models, 
and  in  particular,  allows  us  to  use  realistic  electron  masses  for  FRC  and  magnetic 
reconnection  simulations. 

A  very  simple  multi-dimensional  two-fluid  test  is  a  ring  current  formed  by  strong 
magnetic  field  gradients.  In  this  case  a  magnetic  field  is  embedded  in  a  uniform 
plasma.  The  field  is  positive  within  a  radius  of  0.25  and  negative  otherwise. 
Conducting  walls  are  used  in  the  domain.  The  magnetic  field  gradient  generates 
opposing  electron  and  ion  flows  and  instability  quickly  develops.  We’ve  performed 
this  test  on  both  structured  and  unstructured  grids  to  help  determine  the  grid 
dependence  of  the  solution  for  the  finite  volume  second  order  MUSCL  scheme. 

Figure  3  and  Error!  Reference  source  not  found,  show  solutions  on  an 
unstructured  grid.  The  simulations  were  run  on  8  processors.  The  grid  is 
superimposed  on  the  solution. 

Figure  4  and  Figure  6  show  the  solution  on  a  structured  grid.  The  grid  clearly  plays 
a  role  in  the  final  configuration  of  the  solution,  however,  the  basic  4  fold  symmetry 
is  captured  by  the  solution  on  both  structured  and  unstructured  grids.  It’s  hoped 
that  by  moving  to  a  more  accurate  algorithm  (discontinuous  Galerkin]  the  grid 
dependence  of  the  solution  can  be  mitigated. 

During  the  Phase  I  a  nodal  discontinuous  Galerkin  method  was  implemented  in  Id 
and  2d  in  USim.  Id  nodal  high  order  limiters  were  implemented,  but  not  in  multiple 
dimensions.  The  DG  scheme  works  extremely  well  on  unstructured  grids.  Figure  7 
through  Figure  11  show  a  simple  electromagnetic  pulse  on  an  unstructured  grid 
using  3rd  order  discontinuous  Galerkin  for  Maxwell’s  equations.  The  solution  shows 
great  uniformity  during  expansion  on  the  unstructured  mesh.  We’ve  also 
demonstrated  the  nodal  DG  algorithm  on  a  two-fluid  shock  in  Figure  12  showing 
that  it  can  be  applied  to  multiple  fluids,  Maxwell  equations  and  source  terms 
simultaneously.  In  the  current  version  of  the  code  cell  nodal  values  are  averaged  to 
cell  centered  values  for  plotting  purposes.  In  the  future  nodal  data  will  be  plotted 
node  by  node  so  that  the  solution  can  be  seen  at  the  full  resolution  described  by  the 
DG  approach. 


user:  loverich 
TueJul  3  08:57:18  2012 


Figure  3:  Planar  current  at  time  1  on  an 
unstructured  grid.  Finite  volume  scheme. 
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Figure  4:  Planar  current  at  time  1  on  a 
structured  grid.  Finite  volume  scheme. 
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Figure  5:  Planar  current  at  time  2  on  an 
unstructured  grid.  Finite  volume  scheme. 
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Figure  6:  Planar  current  at  time  2  on  a 
structured  grid.  Finite  volume  scheme. 


Thu  Jan  17  08.05:38  2013 

Figure  7:  Expansion  of  electromagnetic  pulse 
time  1.  3rd  order  nodal  DG  on  unstructured 
mesh. 
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Figure  8:  Expansion  of  electromagnetic  pulse 
at  time  2.  3rd  order  nodal  discontinuous 
Galerkin. 
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Figure  9:  Expansion  of  electromagnetic  pulse 
time  3.  3rd  order  nodal  discontinuous 
Galerkin. 
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Figure  10:  Expansion  of  electromagnetic 
pulse  time  4.  3rd  order  nodal  discontinuous 
Galerkin. 


Figure  11:  Nodal  3rd  order  discontinuous  Galerkin  solution  for  electromagnetic  pulse 
propagation.  Nodal  values  are  averaged  to  cell  center  for  visualization.  Results  are  incredibly 
uniform  despite  the  unstructured  mesh.  No  limiters  were  used  in  this  simulation. 
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Figure  12:  Two-Fluid  solution  using  nodal  discontinuous  Galerkin  method  implemented 
during  the  Phase  I.  During  the  Phase  I  we  implemented  Id  and  2d  nodal  DG  without  multi¬ 
dimensional  limiters.  No  limiters  were  applied  in  this  case  -  however  high  order  nodal 
limiters  are  needed  in  general  and  would  be  a  part  of  the  phase  11  project. 


Figure  13:  A  2D  comparison  of  a  structured  grid  (left:  Ax  =  Ay  =  0.01)  and  an 
unstructured  meshes  (right:  Axmax=  Aymax  =  0.01).  Domain  size  is  1  x  1.  (A,E)  Initial 
conditions  used  for  both  solvers  high  density  =  1.0  (centered  box  0.4  x  0.4) 
superimposed  on  low  density  =  0.1  background).  (B,F)  Close  up  view  of  the  initial 
conditions  for  the  structured  grid  (B)  and  unstructured  mesh  (F).  (C,G)  Solution  to  the 
advection  equation  (density)  at  time  tf  =  0.2  with  flow  speed  of  u  =  v  =  1.  (D,H)  Solution 
to  the  Euler  equations  (density)  at  time  tf=  0.2,  with  zero  initial  velocity  and  constant 
initial  specific  internal  energy  of  e  =  1. 


At  the  University  of  Washington  the  unstructured  finite  volume  method  currently 
under  development  uses  a  k-exact  polynomial  spatial  reconstruction  in  conjunction 
with  a  Runge-Kutta  time  integration  scheme.  This  method  allows  arbitrary  order 


accuracy  for  smooth  solutions  and  drops  to  1st  order  accuracy  in  the  presence  of 
discontinuities.  The  ID  homogeneous  advection.  Burger  and  Euler  equations  are  in 
working  order  along  with  the  2D  advection  equation.  Verification  of  the  2D  Euler 
equations  is  currently  underway.  The  unstructured  mesh  architecture  would  also  be 
ideal  for  use  with  finite  element  and  semi-Lagrange  methods.  Figure  13  shows 
advection  of  a  square  wave  across  both  a  structured  and  unstructured  grid  for 
comparison.  The  solution  shows  good  shape  preservation  despite  the  complex  grid 
structure. 


One  of  the  great  successes  of  the  Phase  I  was  implementation  and  testing  of 
diffusion  based  divergence  cleaning  of  the  electric  field.  Errors  in  charge 


conservation  are  reduced  using  the  equation  —  - c2 V xB  =  —  +  AV  V  •  £  -  — 

dt  £Q  ^  £q  j 

where  the  last  term  is  the  divergence  cleaning  term.  Our  previous  work  with 
electric  field  divergence  cleaning  involved  the  use  of  hyperbolic  divergence  cleaning 
for  the  electric  field  which  worked  very  poorly  in  cases  where  there  was  significant 
charge  separation.  During  this  phase  I  the  parabolic  approach  was  tested  on  a 
simple  problem  (Figure  14)  and  also  tested  on  more  challenging  problems  such  as 
the  GEM  challenge  problems  shown  later  in  this  report.  Parabolic  electric  field 
divergence  cleaning  was  able  to  maintain  small  charge  separation  (small  charge 
conservation  error)  without  destroying  the  solution.  Much  more  testing  of  this 
approach  is  still  needed  and  that  will  be  pursued  during  the  Phase  II  project. 


Figure  14:  Parabolic  cleaning  method  shown  ID  where  the  magnetic  and  electric  fields 
are  initialized  with  random  noise  added  to  their  solution.  Results  show  the  monotonic 
reduction  of  the  divergence  error  of  the  magnetic  and  electric  fields. 


Task  4:  Investigate  techniques  for  solving  chemical  rate 
equations  implicitly 


We’ve  investigated  and  implemented  a  Newton  iteration  combined  with  high  order 
Runge-Kutta  for  generic  stiff  source  terms  including  chemical  reactions  and  multi¬ 
fluid  momentum  and  energy  exchange  terms.  Testing  is  ongoing  and  will  likely  be 
finished  as  part  of  the  Phase  II  of  this  project  or  through  separate  funds. 


Task  5:  Create  test  suite  of  problems  for  verification  and 
validation 

During  the  Phase  I  we  pursued  verification  and  validation  tests  for  the  new  semi- 
implicit  time  integration  scheme.  Using  this  scheme  we  were  able  to  obtain  results 
using  realistic  electron  to  ion  mass  ratio  while  maintaining  high  temporal  accuracy 
for  FRC  formation  and  magnetic  reconnection. 

One  of  the  more  challenging  problems  is  FRC  formation  in  axisymmetric  geometry. 
Initial  conditions  for  this  simulation  were  provided  by  Richard  Milroy  at  the 
University  of  Washington.  Table  1  through  Table  4  show  results  of  5  moment  two- 
fluid  FRC  formation  using  realistic  electron  to  ion  mass  ratio.  In  order  to  get  the 
simulations  to  complete  on  time  we  dropped  the  speed  of  light  by  a  factor  of  100, 
though  running  these  simulations  with  the  exact  speed  of  light  is  quite  possible.  In 
order  to  run  the  simulation  we  first  approximated  coil  fields  by  running  our 
Maxwell's  equation  solver  to  a  quasi-steady  state  (no  interaction  between  fluids  and 
fields),  at  which  point  the  full  two-fluid  simulation  began.  The  FRC  chamber  is 
initially  filled  with  an  ionized  gas  at  2ev  with  a  sinusoidal  field  applied  at  the 
boundary  and  crowbar  at  40  microseconds.  Table  1  shows  the  solution  near  the 
point  where  the  EM  field  reaches  steady  state,  just  as  the  two-fluid  simulation  is 
turned  on.  Table  2  shows  the  solution  as  the  field  reverses  on  the  boundary  and  the 
plasma  is  pushed  towards  the  axis  and  begins  to  collapse  towards  Z=0.  Table  3 
shows  the  solution  at  28  microseconds  while  the  FRC  is  still  compressing  and 
collapsing.  Table  4  shows  the  solution  at  crowbar  and  the  FRC  has  reached  a  relative 
steady  state.  We’ve  attempted  these  same  simulations  with  10-moment  ions,  but  no 
success  so  far.  We  may  still  need  to  work  out  issues  with  the  axisymmetric  source 
terms  -  however,  we  believe  simulations  in  3D  (instead  of  axisymmetric  geometry) 
with  10  moment  ions  would  probably  work. 


Table  1:  FRC  at  6.4  microseconds 

Z  magnetic  field 


Pseudocolor 
Var:  flulds/em  5 
-0.6446 


(0.4376 
0.2305 

0.02348 

-0.1836 
Max:  0.6446 
Min: -0.1836 


J? 


Magnetic  field  lines 


Streamline 
Var:  Speed 

0.1843 


0.1432 


0.1022 


0.06119 


0.02017 
Max:  0.1843 
Min:  0.02017 


N 


J? 


Electron  number  density 


R 


Table  2:  FRC  at  16  microseconds 

Z  magnetic  field 


Pseudocolor 
Var:  flulds/em  5 
—  0.6444 


Max:  0.6444 
Min: -0.21 31 


3.0 


1.0 


N  0.0 


-1.0 


-2.0 


-3.0 


R 


Magnetic  field  lines 


Electron  number  density 


R 


Table  3:  FRC  at  28  microseconds 

Z  magnetic  field 


Pseudocolor 
Var:  flulds/em  5 
-6.8534 


Max:  0.8534 
Min: -0.3334 


R 


Magnetic  field  lines 


Electron  number  density 


Pseudocolor  ? 
Var:  ne 

-5.846e+21 


I 


4.385e+21 
2.925e+2l2 . 0 
- 1 ,464e+21 


^^™-3.322e+18 
Max:  5.846e+21  J  .  0 
Min:  3.322e+18 


A  very  common  problem  for  testing  two-fluid  effects  is  the  GEM  challenge  magnetic 
reconnection  problem.  These  simulations  were  performed  with  both  5  moment 
two-fluid  and  10  moment  ions  with  5  moment  electrons.  There  is  one  key 
distinction  between  the  two  models  during  magnetic  reconnection,  i.e.  the  transfer 
of  kinetic  energy  to  thermal  energy.  With  10  moment  the  kinetic  energy  of 
reconnection  is  preferentially  transferred  to  Pxx,  the  ion  pressure  in  the  X  direction. 
Using  5  moment  ions  it’s  assumed  that  Pxx,  Pyy  and  Pzz  are  equal  so  the  energy  is 
also  transferred  into  the  Y  and  Z  components  of  the  pressure  tensor.  The  result  is 
that  in  the  10  moment  simulations  we  see  less  expansion  of  the  plasmoids  in  the  y 
direction  as  Pyy  remains  smaller  than  in  the  5  moment  ion  case.  These  simulations 
were  performed  using  the  new  semi-implicit  scheme  with  realistic  electron  to  ion 
mass  ratio  and  serve  as  a  verification  test  of  our  algorithm  for  both  the  5  moment 
and  10  moment  models  in  2D.  Figure  15,  Figure  16,  Figure  17  show  the  5  moment 
solution  at  various  times.  Figure  18  and  Figure  19  show  the  10  moment  solutions  as 
various  times.  Figure  20  shows  the  reconnected  flux  for  both  the  10  moment  and  5 
moment  solutions.  In  the  5  moment  solution  magnetic  island  formation  is  observed, 
the  island  then  merges  with  the  main  island.  After  island  merger  reconnection  rates 
from  both  the  10  moment  and  5  moment  solutions  agree. 


Table  4:  Frc  at  40  microseconds 
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Figure  15:  Electron  Z-momentum  for  5  moment  two-fluid  reconnection  at  time  0. 


Figure  16:  Electron  Z  momentum  for  5  moment  two-fluid  reconnection  at  30  ion  cyclotron 
periods  using  realistic  electron  to  ion  mass  ratio. 
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Figure  17:  Ion  density  for  5  moment  two-fluid  reconnection  at  20  ion  cyclotron  periods  using 
realistic  electron  to  ion  mass  ratio. 


Figure  18:  Electron  Z  momentum  for  10  moment  two-fluid  reconnection  at  time  30  ion 
cyclotron  times  using  realistic  electron  to  ion  mass  ratio. 


Figure  19:  Ion  density  for  10  moment  two-fluid  reconnection  at  time  20  ion  cyclotron  periods 
using  realistic  electron  to  ion  mass  ratio. 


Reconnected  flux  vs  time  for  two-fluid  models 


Figure  20:  Reconnected  flux  using  5  moment  and  10  moment  models  with  realistic  electron  to 
ion  mass  ratio. 


